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Abstract 


In this study, we present an improved physical model to retrieve snowfall rate over 
land using brightness temperature observations from the National Oceanic and 
Atmospheric Administration’s (NOAA) Advanced Microwave Sounder Unit-B (AMSU- 
B) at 89 GHz, 150 GHz, 183.3±1 GHz, 183.3±3 GHz, and 183.3±7 GHz. The retrieval 
model is applied to the New England blizzard of March 5, 2001 which deposited about 75 
cm of snow over much of Vermont, New Hampshire, and northern New York. 

In this improved physical model, prior retrieval assumptions about snowflake shape, 
particle size distributions, environmental conditions, and optimization methodology have 
been updated. Here, single scattering parameters for snow particles are calculated with 
the Discrete-Dipole Approximation (DDA) method instead of assuming spherical shapes. 
Five different snow particle models (hexagonal columns, hexagonal plates, and three 
different kinds of aggregates) are considered. Snow particle size distributions are 
assumed to vary with air temperature and to follow aircraft measurements described by 
previous studies. 

Brightness temperatures at AMSU-B frequencies for the New England blizzard are 
calculated using these DDA calculated single scattering parameters and particle size 
distributions. The vertical profiles of pressure, temperature, relative humidity and 
hydrometeors are provided by MM5 model simulations. These profiles are treated as the 
a priori data base in the Bayesian retrieval algorithm. In algorithm applications to the 
blizzard data, calculated brightness temperatures associated with selected database 
profiles agree with AMSU-B observations to within about ±5 K at all five frequencies. 
Retrieved snowfall rates compare favorably with the near-concurrent National Weather 
Service (NWS) radar reflectivity measurements. The relationships between the NWS 
radar measured reflectivities Z e and retrieved snowfall rate R for a given snow particle 
model are derived by a histogram matching technique. All of these Z e -R relationships fall 
in the range of previously established Z e -R relationships for snowfall. This suggests that 
the current physical model developed in this study can reliably estimate the snowfall rate 
over land using the AMSU-B measured brightness temperatures. 
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1. INTRODUCTION 


Falling snow is an important component of the global water cycle. Heavy snowfall 
can disrupt transportation flow and cause subsequent severe flooding. Snow packs 
accumulated on the ground serve as a reservoir of water for agriculture and hydroelectric 
power generation. Moreover, falling snow that persists as snow cover over land can affect 
earth energy balance through a change in the surface albedo. Understanding extratropical 
precipitation is critical for improving the prediction capability of regional and large-scale 
climate models for the water cycle. For the last several decades, ground-based radars and 
snow gauges have been used to monitor snowfall rate. However, spatial coverage of radar 
and snow gauge networks outside of the USA, Europe, and Japan is sparse. Snowfall rate 
measurements from space can overcome this spatial sampling limit and provide data sets 
necessary for the improvement of weather forecasting, hydrological and climate research. 

While satellite-based rain rate estimates are reliable and operational (Olson et al. 
1996; Ferraro el al. 2005; Olson et al. 2006), the measurement of snowfall rates from 
space is a relatively new field (Chen and Staelin 2002; Kongoli et al 2003; Skofronick- 
Jackson et al. 2004; Liu 2004, Noh et al. 2006). There are two major challenges 
associated with retrieving snowfall rates (a) adequately representing and retrieving the 
complex macro and microphysical features of snow clouds, and (b) distinguishing surface 
features from atmospheric signatures. The measurement of snowfall within the 
atmosphere has been difficult using radiometers that operate at frequencies below 100 
GHz where the atmosphere is relatively transparent and the surface emissivity can 
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produce brightness temperatures (Tb) expected from precipitating clouds. Indeed, 
retrievals of snow pack rely on the 19 and 37 GHz channels (Foster et al. 2005). 
Recently, falling snow retrievals have been derived from spaceborne microwave 
radiometry over oceanic regions where the measurements are not affected by surface 
snow (Weinman and Hakkarinen 1990; Liu and Curry 1996; Schols et al. 1999; Bennartz 
and Petty 2001). 

Passive microwave radiometers operating at frequencies near 183 GHz, such as 
NOAA’s Advanced Microwave Sounder Unit-B (AMSU-B), and those on proposed 
satellites like the Global Precipitation Mission (GPM) and the Meteorological 
Operational Mission (METOP) polar orbiting satellites, can minimize the surface 
emission problem over land because water vapor absorption effectively masks the surface 
emission. For example, AMSU-B operating at 89 GHz, 150 GHz, 183. 3±7, ±3, and ±1 
GHz has been employed to estimate frozen hydrometeors (Chen and Staelin 2002; 
Kongoli et al 2003; Skofronick-Jackson et al. 2004; Noh et al. 2006). Kongoli et al. 
(2003) derived snowfall over land from the AMSU-B Tbs using empirical relationships. 
Chen and Staelin (2003) employed a neural network technique in their snowfall 
retrievals. Although such empirical relationships or statistical techniques may be 
operationally useful, physical models are needed to understand how the measured Tbs 
depend on various surface and atmospheric parameters. 

Skofronick-Jackson et al. (2004) developed a physical model at millimeter-wave 
frequencies from which snowfall rates over land could be inferred. While their snowfall 
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rate retrievals were qualitatively validated, the retrievals are greatly enhanced by 
improving the assumptions. 


• That study approximated the single scattering properties of snow crystal as equal 
volume (V) to area (A) ratio spheres (Grenfell and Warren 1999). In the IR/UV 
frequency regions, the methods using the equal-V/A spheres generally do not 
provide accurate estimates of the asymmetry factors (Takano and Liou 1988; 
Grenfell and Warren 1999). Recent analyses by Liu (2004) and Kim (2004) 
demonstrated the limitation of such an approach at microwave frequencies. A 
refined technique to calculate the scattering properties of non-spherical snow 
particles is employed herein. 

• In addition, the particle size distribution (PSD) of equal-V/A spheres employed in 
the early study was simplified over the vertical height of the cloud. In this work, 
PSDs are taken from in-situ observations of snowflakes over the vertical domain 
of clouds. 

• Another simplifying assumption in the previous study was the number of 
atmospheric and hydrometeor profiles used in the retrievals. Only 1 1 and 36 
representative profiles of relative humidity and snow water content, respectively, 
were extracted from mesoscale model simulations. Additionally, the model 
simulation resolution from the earlier study was 40 km x 40 km, which is much 
larger than the AMSU-B footprints (~ 16 km x 16 km at nadir). To obtain more 
realistic retrievals, the number and diversity of atmospheric profiles needed to be 
expanded and the resolution of the simulated data improved. 
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• Further, the method of Skofronick-Jackson et al. (2004) was to find a final profile 
that minimized the differences between the observed Tbs and calculated Tbs 
without any consideration of the correlation between model variables (e.g. very 
low relative humidity was retrieved with heavy snowfall). The current study uses 
a Bayesian method for optimizing the retrievals and correlations between 
variables are considered in retrievals. Noh et al.(2006) also developed a Bayesian 
method based algorithm to retrieve the froze precipitation retrievals using high 
frequency microwave satellite data. However, their study did not consider error 
correlations between different channels so that only diagonal terms of the error 
covariance matrix were employed in the Bayesian algorithm. In this study, the 
correlations between the retrieved variables are considered. Modeling errors are 
obtained from systematic Tb sensitivity tests, based upon the uncertainties of 
several parameters, including particle size, water vapor content, ice content, and 
surface emissivity. 

The purpose of this study is to improve the retrieval algorithm for estimating 
snowfall rate with millimeter- wave channels. The major improvements from this study 
are given below. First, the Discrete-Dipole Approximation (DDA) method (Purcell and 
Pennypacker 1973; Draine and Flatau 2003) is employed to calculate more precisely the 
single scattering parameters of nonspherical snow crystals in radiative transfer 
calculations. Second, snow particle size distributions are based on in-situ observations by 
Houze et al (1979), Lo and Passarelli (1982), and these distributions are allowed to vary 
with height. Third, an a priori data base containing appropriate and realistic simulated 
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profiles is created for snowfall parameter retrievals within a Bayesian methodology. In 
addition, an error analysis and sensitivity tests are perfonned. 

The paper is organized as follows: The blizzard snow storm case considered in this 
study is briefly described in section 2; radar and AMSU-B observations as well as MM5 
simulated cloud profiles are also presented. Details of the a priori database are explained 
in Section 3. The retrieval methodology including the Bayesian method is described in 
section 4. Retrieval results and uncertainty analysis are presented in Section 5. Finally, 
summary and conclusions are given in Section 6. In the Appendix, the sensitivity of the 
AMSU-B radiometer to humidity and precipitation profiles and surface snow coverage 
are examined by calculating Jacobians (after Bauer and Mugnai 2003). 

2. SNOWSTORM CASE: Northeastern Blizzard on March 5, 2001 

The snowstorm considered in this study is the blizzard of March 5-6, 2001 over the 
Northeastern United States. This is the same blizzard case reported by Skofronick- 
Jackson et al. (2004). The blizzard deposited about 75 cm of snow over much of Vermont, 
New Hampshire, and northern New York during its lifetime. Figure 1(a) shows a 
composite of the National Weather Service (NWS) operational weather radar reflectivity 
March 5, 2001 at 2300 UTC. The NWS radar reflectivity data were averaged over 16 km 
x 16 km grid template to match the finest spatial resolution of the AMSU-B channels. 
The maximum reflectivity in the smoothed radar reflectivity data over land is ~40 dBZ. 
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Shown in Figures 1 (b)-(f) are the Tbs measured by the AMSU-B on the National 
Oceanic and Atmospheric Administration (NOAA)-15 satellite at 2302 UTC on March 5, 
2001. The AMSU-B is a cross-track scanner with an angular swath of ±50°. The 
southeast to northwest intense snow band is at about 35° off nadir. It is noted that the 89 
GHz channel shows ambiguity in distinguishing snow in the atmosphere from other 
surface features on the ground such as lakes, rivers, and snow on the ground. On the other 
hand, the 183 GHz channels (especially those closest to the water vapor line center) can 
screen the surface effect except in the driest atmospheric conditions such as those found 
in the arctic. 


3. GENERATION OF A PRORI DATA BASE 
3. 1 Mesoscale model simulation 

In the retrieval algorithm, vertical profiles of pressure, temperature, relative 
humidity, cloud ice, cloud water, and hydrometeor profiles are derived from the MM5 
model (Grell et al. 1994) which was used to the Northeastern Blizzard at 4 km resolution. 
The MM5 model was initialized at 00 UTC on March 5, 2001 and the model integration 
was perfonned for a period of 24 hours. The simulation domain with 36 km resolution 
was nested to the domain with 12 km resolution and then to the domain with 4 km 
resolution. The model domain was centered at 35° N and 70° W and consisted of 231 x 
231 grid points at 4 km resolution. These profiles were linearly averaged to 16 km x 16 
km resolution to match the AMSU-B foot print size at nadir. 
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The Goddard scheme (Tao and Simpson 1993) was used for ice microphysics 
parameterizations in the simulations. All the temperatures near the surface over New 
England were below -2 °C and the observations are at 19:00 local time; therefore, the 
precipitation was either snow or graupel over land. Graupel generated by the MM5 model 
over land was converted to snow of equivalent water content in this study because the 
National Lightning Detection Network (NLDN) observed some lightning off-shore, but 
not over New England. 

In addition to the MM5-generated profiles at 16 km x 1 6 km resolution, assumed 
fractional surface snow cover values of 0, 0.1, 0.2, ..., 1.0 are used as input to forward 
radiative transfer calculations to adjust the surface emissivity. Thus, the surface varies 
from 100 % mixed bare soil/frozen soil/winter forest to 100 % snow cover. Employing 
these input profiles and calculated Tbs in an a priori database, a Bayesian algorithm is 
used to retrieve atmospheric environment and hydrometeor profiles and surface snow 
cover fraction using AMSU-B observed Tbs. Details of the scattering parameterizations, 
particle size distributions, snow particles’ orientation, radiative transfer calculations, and 
the Bayesian retrieval framework are presented in the following sections. 

3.2 Single scattering parameter of snow particles 

The Discrete Dipole Approximation (DDA) method was used to compute single 
scattering parameters of various idealized nonspherical snow crystals which were then 
used in the Tb calculations. The DDA method is a flexible technique for calculating the 
electromagnetic scattering and absorption by particles with arbitrary shapes and 
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composition (Draine 1988; Mishche nk o et al. 2000). The DDA treats the actual particle 
as an array of dipoles. Each of the dipoles is subject to an electric field which is the sum 
of the incident wave and the electric fields due to all of the other dipoles. Through the 
solution of the electric field at each dipole position, the scattering and absorption 
properties of the particle are obtained. The DDA replaces a solid particle with an array of 
point dipoles occupying positions on a cubic lattice, and the lattice spacing must be small 
compared to the wavelength of the incident radiation. Therefore, the DDA method 
requires large computer storage and CPU time. The technique is not well suited for 
particles with very large complex refractive indices because it requires much narrower 
distance between dipoles, thus requiring much larger memory size. This study employs 
the DDA codes developed by Draine and Flatau (2003). 

The five idealized snow crystal models considered in this study are shown in Figure 
2: hexagonal columns (HC), three types of snow aggregates composed of two cylinders 
(C2), three cylinders (C3) and four cylinders (C4), and hexagonal plates (HP). 
Singlescattering parameters of HCs, C2s, C3s, C4s, and HPs were previously calculated 
by Kim (2006) for the aspect ratio (the ratio of small dimension to large dimension of ice 
crystals) ~ 0.1. 

For this work we changed the aspect ratio to one that had a thinner column based on 
the observations in the previous studies (Auer and Veal 1970; Heymsfield 1972). 
Following Auer and Veal (1970) and Heymsfield (1972), the relationship between 
diameter (D) and length (F) for hexagonal columns is given by 
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D = 0.197 xL 0414 [mm] (1) 

The aggregates (C2s, C3s, and C4s) are modeled with two, three, and four circular 
cylinders having the same aspect ratio as a hexagonal column. 

The thickness (T) and width (W) relationship for hexagonal plates follows Auer and 
Veal (1970) and is given by 

T = 0.048 xW“ [mm] (2) 

The density of ice for each cylinder and plate is assumed to be the same as pure ice (0.91 
g/nr ). Dielectric constants are calculated with formulas given by Matzler and Wegmuller 
(1987). 

3.3 Orientation of snowflakes 

Single scattering parameters of nonspherical particles strongly depend on the 
orientation of the particles in the snowfall. Vivekanandan et al. (1994, 1999) showed that 
differential reflectivities of snow particles measured by polarimetric radar observations 
were close to zero. Differential reflectivity is the ratio of the horizontal copolar return to 
the vertical copolar return and can be interpreted as the reflectivity weighted mean-axis 
ratio of the precipitation particle in the radar resolution volume. Thus it is a good 
indicator of orientation of particles. This implies that snow particles are randomly 
oriented in dry snow regions of various precipitating systems. For this reason, snow 
particles are assumed to be randomly oriented in our analysis. 
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In addition, the atmosphere was known to be strongly turbulent in the snowstorm 
under study. Shown in Figure 3 are the profdes of wind speed and Richardson number 
observed by NWS instrumentation in Brookhaven (41° N, 73° W), New York at 00 UTC 
on 6 March 2001, an hour after the NOAA-15 AMSU-B observations; This figure shows 
that wind speed was greater than 10 m/sec with strong wind shear at most altitudes. 


The Richardson number, Ri (Bluestein 1993), which is a measure of the importance 
of buoyancy forces to inertial accelerations, is defined as follows: 



( 3 ) 


where 9 is the potential temperature and v is the wind speed. When Ri is small 
(<0.25), the flow becomes turbulent (Bluestein (1993)). Small Ri numbers were observed 
during this storm, suggesting that the atmosphere below 10 km height was turbulent and 
the particles are randomly oriented. 


It should be stressed that this case study is very limited and the random orientation of 
snow particles assumed here may not be generalized to all snowstonns. For example, 
Hogan et al. (2002) analyzed simultaneous aircraft and polarimetric radar data for a 
warm- frontal mixed-phased cloud. By analyzing the differential reflectivity given by 

f Z ' 

Z dr = 101og 10 — — dB , where Zh and Zy are reflectivity factors measured at horizontal 

V Z v ) 

and vertical polarizations, respectively, they showed the possibility of horizontal 
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alignment of the ice crystals in the region of embedded convection where high 
concentrations of small crystals were observed. More comprehensive observations of the 
orientation of snow particles using radar and passive microwave radiometer 
measurements, in-situ snow crystal microphysics samples, atmospheric stability 
observations, and wind profiles from high resolution soundings are necessary to clarify 
the particle orientation issue, and to gain a more complete understanding of falling snow 
retrievals. 

3.4 Snow particle size distributions 

In-situ microphysics data describing snow crystal shapes and size distributions were 
not available during this snow storm. Therefore, we employ snow particle size 
distributions (PSDs) measured by Houze et al. (1979) and Lo and Passarelli (1982) for 
midlatitude winter storms. 

According to these studies, the snow PSDs are represented by exponental functions 
given by 

N(D) = N 0 exp(-AD) [mm 1 in' 3 ] (4) 

where D is the large dimension of snow particle. 

By fitting a curve to the data shown in Houze et al. (1979), A is assumed to follow: 

T 

A = lO~ Tl [mm' 1 ] (5) 

where T(°C) is the air temperature. It is noted that Eqn.(5) can be used to describe the X 
versus air temperature relationship shown in Lo and Passarelli (1982). 

3.5 Falling velocity of snow particles 
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To convert snow water contents to snowfall rates, the following equation by 
Rutledge and Hobbs (1983) for terminal velocity (V t ) of snow particles is employed: 

V t (D) = 1.139 x D° 11 x (^-) 0 ' 4 [m/sec] (6) 

where P is atmospheric pressure for a given altitude and Po is a reference value equal to 

p 

1000 mb (Rutledge and Hobbs 1983). The factor (^-) 04 allows for the change in fall 

speed with air pressure (Foote and DuToit, 1969), such that as the pressure decreases 
with altitude, the particles fall faster. 

Then snowfall rates are calculated using 

RR(mm/hr) = 6 ttx 0.91 x 10 -4 jD^V t (D)N(D)dD (7) 

where D e ff [mm] is the diameter of an equal-mass ice sphere corresponding to a given 
snow particle. The factor of 0.91 is multiplied in Eqn.(7) to convert ice snow to an 
equivalent volume of liquid water. 

3.6 Radiative transfer calculations 

The radiative transfer model used to compute Tbs from a hydrometeor profile is the 
delta-Eddington model (Weinman and Davis 1978, Thomas and Stamnes 1999, Kim et al. 
2004). Vertical profiles of pressure, temperature, humidity and hydrometeors are 
provided by the MM5 model simulations described in Section 3.1. Single scattering 
coefficients, asymmetry factors, and single scattering albedos of snow particles are 
calculated with scattering models described in Section 3.2, based upon the particle size 
distributions described in Section 3.4. Gaseous absorption coefficients for oxygen and 
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water vapor are obtained from the millimeter propagation model (MPM) of Liebe et al. 
(1992). 

The radiative transfer calculations also require a knowledge of the emissivity of the 
variable surface features, including snow cover. The boundary conditions were 
determined partially by the accumulated antecedent snow for which the emissivities e s 

for deep dry snow at a 35° viewing angle are 0.64, 0.724, and 0.8 at 89, 150, and 183 
GHz, respectively (Hewison and English 1999). The emissivity used in the radiative 
transfer model is a weighted mean of the emissivity of snow cover s s and £ 0 , where s 0 
is an average of the emissivities of bare soil, frozen soil, and winter forest/conifer 
surfaces. The value of s 0 is 0.98 for all frequencies based on the observations shown in 
Hewison et al. (2001). The effective emissivity is thus 

S = f£ s +(l-f)e 0 (8) 

where f is the snow cover fraction that is assigned 11 discrete values: 0.0, 0.1, 0.2, 0.3, 
0.4, ...., 1.0. 


4. RETRIEVAL METHODOLOGY 


A Bayesian inversion technique is used to retrieve falling snow profiles. An a priori 
data base was built using the MM5 generated atmospheric and hydrometeor profiles. For 
a given snowfall profile, different snow particle shapes (5 kinds) and various surface 
emissivity values (11 snow cover fractions) were considered in the radiative transfer 
calculations. That is, 1 1 different Tbs were generated at each frequency for a given 
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snowfall profile and a given snow particle type. In-situ microphysics observations of 
snow particle habits were not available; therefore, the five snow particle models 
discussed in Section 3.2 are used in the Tb calculations. However, for a given retrieval 
application to the AMSU-B data, only one shape was allowed at a time due a lack of 
information about the number and types of particle shapes for this blizzard case. 
Retrievals assuming different particle shapes are compared. 

Following Olson et al. (1996) and Moreau et al. (2003), if it is assumed that the 
errors in the observations and the simulated observations are Gaussian and uncorrelated 
then the “best” estimate of state vector x, given the set of observations y 0 , is the 
expected value 

x a = J J J xexp {-°* 5 [y 0 - y s ( x )l T (° + s ) " I [y 0 -y s ( x )]} p a( x = x true) /A dx ( 9 ) 

where A is a normalization factor, 

A = n - j ex P^ _0 - 5 [yo-y s ] T (0 + S) '[yo-y s (x)]}P a (x = x true )dx. (10) 

Here, y s (x) are simulated Tbs and P a is the a priori probability that x is the true state 
profile of the atmosphere. The O and S are the observed and simulated Tb error 
covariance matrices, respectively. 

The inherent uncertainty is given by the integration of departures of the best estimate 
state vector x a from those contained in the database: 

E (x) = J J... J(x -x a ) 2 e xp{-°.5[y 0 - y s (x)] T (O + S) ~'[y 0 -y s (x)]}P a (x = x true )/A dx 

(ID 
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When generating the a priori database, the evaluation of the associated modeling 
uncertainties is difficult but important for both retrievals and error estimation. Tb 
uncertainties were included in the error covariance matrix representing the intrin s ic 
variability of the generated a priori database. Tassa et al. (2006) proposed a methodology 
for taking into account cloud-radiation database related modeling uncertainties for 
precipitation retrieval from the TRMM Microwave Imager (TMI). In their study, the 
modeling errors were obtained from systematic Tb sensitivity tests, based upon the 
assumed uncertainties in PSD, atmospheric temperature, ice content, sea surface wind 
speed, viewing angle, melting phase, and particle shape. 

Following Tassa et al. (2006), modeling uncertainties associated with the generation 
of a priori database were calculated in this study for snowfall retrievals using AMSU-B 
measurements. The uncertainties were computed through sensitivity analyses aimed at 
evaluating the impact of various bulk cloud/radiative parameters on the resulting 
simulated radiance distributions. We assume there are four fundamental sources of error: 
errors in calculating surface emissivity, errors in calculating scattering parameters, errors 
in water vapor profiles, and errors in radiative transfer. We address these sources of error 
as follows: 

(a) Regarding surface emissivity, we add Gaussian-distributed errors to our 
computed emissivity for each emissivity calculation we make in generating the database. 
The standard deviation of emissivity errors assumed in this study is 0.05. 
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(b) Even though we consider the Bayesian estimates for each particle habit 
separately, there still could be errors in each calculation of the single scattering 
parameters and the size spectra of the hydrometeors that can have a very large impact on 
snow radiative properties (Panegrossi et al. 1998, Viltard et al. 2000). In this study, we 
added Gaussian-distributed errors to the mass median diameter of the PSDs described by 
Eqn.(4). The standard deviation for the percentage error was assumed to be 50 %. 

(c) Since AMSU-B brightness temperatures are sensitive to water vapor profiles, we 
added Gaussian-distributed errors to the MM5 generated water vapor content profiles. 
The standard deviation of the percentage error was assumed to be 10 % and the error was 
assumed to be same for all levels. 

(d) Because we did not use rigorous a Monte Carlo radiative transfer model because 
of increased computation time required for these calculations, there will be some 
uncertainty in radiance calculations that increases in proportion to the effective scattering 
in the vertical column. Therefore, we calculate a “cloud -free” background Tb and then 
calculated the difference between the background and cloudy atmosphere Tb. This Tb 
difference represents the effective scattering by snow in the vertical column. We added a 
Gaussian-distributed percentage error to the cloudy atmosphere Tb in proportion to a Tb 
difference. The standard deviation of the percentage error was assumed as 8 %. For 
example, if the original Tb was 100 K below the background value, then typically we 
would add or subtract 8 K to the Tb. Since the Tb depressions at 89 GHz and 150 GHz 
relative to the background can be affected by the low surface emissivity of snow on the 
ground, errors from source (d) were calculated by assuming surface snow fraction was 
zero and isolating the Tb depression due to scattering in the atmospheric column. 
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Errors only from (a) first, and then (b) only, and then (c) only, and then (d) only, and 
then all four error sources together were calculated. The error covariances due to errors 
from (a), (b), (c), and (d) are shown in Tables 1-4, respectively. The error covariances 
calculated with (a), (b), and (c) error sources considered together are shown in Table 5. 
Finally, each component of Table 4 and Table 5 are summed for the total error covariance 
(Table 6) which was applied in Eqn. (9) for snow retrievals. 


Noh et al. (2006) neglected the error correlations between different channels by 
assuming the off-diagonal tenns of error covariance matrix to be zero. However, this 
study demonstrates that the modeling error correlations are significantly large between 
different channels of AMSU-B and emphasizes that they should not be neglected in 
precipitation retrievals. 


2 

Table 1. Computed error covariances (K ) due to error source (a) surface emissivity 



89 GHz 

150 GHz 

183.3±7 GHz 

183. 3±3 GHz 

183. 3±1 GHz 

89 GHz 

39.32 

16.82 

0.33 

0.00 

0.00 

150 GHz 

16.82 

7.78 

0.17 

0.00 

0.00 

183.3±7GHz 

0.33 

0.17 

0.005 

0.00 

0.00 

183.3±3GHz 

0.00 

0.00 

0.00 

0.00 

0.00 

183.3±lGHz 

0.00 

0.00 

0.00 

0.00 

0.00 
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Table 2. Computed error covariances (K“) due to error source (b) snow particle size 



89 GHz 

150 GHz 

183.3±7 GHz 

183. 3±3 GHz 

183. 3±1 GHz 

89 GHz 

19.45 

38.43 

6.74 

-2.37 

-1.15 

150 GHz 

38.43 

89.71 

23.57 

-0.28 

-1.09 

183.3±7GHz 

6.74 

23.57 

15.58 

6.59 

1.67 

183.3±3GHz 

-2.37 

-0.29 

1.77 

5.57 

6.59 

183.3±lGHz 

-1.15 

-1.09 

1.67 

1.77 

0.61 
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Table 3. Computed error covariances (K“) due to error source (c) water vapor amount 



89 GHz 

150 GHz 

183.3±7 GHz 

183. 3±3 GHz 

183. 3±1 GHz 

89 GHz 

8.82 

11.72 

0.92 

-3.16 

-4.71 

150 GHz 

11.72 

16.42 

1.82 

-4.31 

-6.60 

183.3±7GHz 

-0.92 

1.82 

1.29 

0.09 

-0.41 

183.3±3GHz 

-3.16 

-4.31 

0.09 

1.58 

2.07 

183.3±lGHz 

-4.71 

-6.60 

-0.41 

2.07 

3.00 
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Table 4. Computed error covariances (K“) due to error source (d) radiance computation 



89 GHz 

150 GHz 

183.3±7 GHz 

183. 3±3 GHz 

183. 3±1 GHz 

89 GHz 

0.86 

1.62 

0.84 

0.57 

0.89 

150 GHz 

1.62 

3.15 

1.61 

1.04 

1.56 

183.3±7GHz 

0.84 

1.61 

0.88 

0.59 

0.89 

183.3±3GHz 

0.57 

1.04 

0.59 

0.47 

0.75 

183.3±lGHz 

0.89 

1.56 

0.89 

0.75 

1.25 
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Table 5. Computed error covariances (K“) due to error sources (a), (b), and (c) 



89 GHz 

150 GHz 

183.3±7 GHz 

183. 3±3 GHz 

183. 3±1 GHz 

89 GHz 

70.87 

66.79 

2.62 

-7.92 

-6.29 

150 GHz 

66.79 

98.68 

9.99 

-10.22 

-7.66 

183.3±7GHz 

2.62 

9.99 

5.77 

3.23 

1.78 

183.3±3GHz 

-7.92 

-10.22 

3.23 

5.98 

3.88 

183.3±lGHz 

-6.29 

-7.66 

1.78 

3.88 

3.54 


2 

Table 6. The error covariances (K ) employed in the retrievals in this study 



89 GHz 

150 GHz 

183.3±7 GHz 

183. 3±3 GHz 

183. 3±1 GHz 

89 GHz 

71.73 

68.41 

3.46 

-7.35 

-5.4 

150 GHz 

68.41 

101.83 

11.60 

-9.18 

-6.1 

183.3±7GHz 

3.46 

11.60 

6.57 

3.82 

2.67 

183.3±3GHz 

-7.35 

-9.18 

3.82 

6.45 

4.63 

183.3±lGHz 

-5.4 

-6.1 

2.67 

4.63 

4.79 


Table 6 shows that modeling errors are strongly correlated between 89 GHz and 150 
GHz channels and between 183.3±3 and 183.3±1 GHz channels and correlations 
coefficients are 0.80 and 0.83, respectively. It is noted that water vapor uncertainties 
result in Tb errors that are significantly anti-correlated at the lower and higher channel 
frequencies. That is, more water vapor increases the Tb at 89 GHz but decreases the Tb at 
183. 3±1 GHz. Table 5 shows that the largest Tb error covariances are associated with the 
89 and 150 GHz channels, primarily due to assumed uncertainties in modeled surface 
emissivities and PSDs. 
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5. RETRIEVAL RESULTS 


5. 1 Falling Snow Retrievals 

Using the algorithm previously described, retrievals were performed for the March 
2001 blizzard case shown in Figure 1. Each retrieved profile contains vertical 
distributions of temperature, relative humidity, and snow water content. The lowest 
altitude snow water content (at 20.0 m) is converted to a melted precipitation rate using 
the procedure described in Section 3.5 (Falling Velocity of Snow Particles), as shown in 
Figure 4(a). Since the retrieval method was designed for applications over land, the 
oceanic regions are masked. The snowflake shape employed in the retrievals shown in 
this figure is HC (Fig. 2). Retrieval results from different snow particle models are 
compared later in this section. As seen in Figure 4(a), the spatial distribution of retrieved 
snowfall rates is similar to the radar reflectivity observations displayed in Figure 1(a). 
Retrieved water equivalent snowfall rates reach a maximum of ~ 4.3 mm/hr, which helps 
to explain the heavy snowfall accumulations (up to ~75 cm) reported during this storm. 

In order to understand the retrievals on a Tb basis, the spatial distributions of 
differences between the Bayesian retrieved Tbs and the AMSU-B observations at all five 
frequencies of the stonn region are shown in Figure 5. Results show that the computed 
Tbs agree with the observations within ±5 K over most of the stonn region at all of the 
AMSU-B frequencies. 

Using Eqn.(ll), the inherent uncertainty of retrievals is calculated by integrating 
departures of profiles contained in the a priori database from the retrieved profile. The 
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distribution of error standard deviations related to the retrieved snowfall rates is shown in 


Figure 4(b). It may be inferred that the uncertainty of the retrieved snowfall rates in this 
study range between 0.01 and 2 mm/hr at 0.02 km altitude. As shown in Figure 4(c), the 
error standard deviations increase with the retrieved snowfall rates. 

Similar to Figure 4, the retrievals are repeated for each of the five snow particle 
shapes shown in Figure 2. To compare retrieved snowfall rates using the different snow 
particle models, a pixel area matching technique similar to that described in Calheiros 
and Zawadski (1987) is used to relate the radar reflectivity Z e (mm m' ) observed by 
NWS radars over land to the retrieved snowfall rate, R at 0.02 km altitude above the 
surface. The procedure selects a number of pixels that exceed a given snowfall rate and 
the same number of pixels that exceeded a particular radar reflectivity. Threshold values 
of each of these quantities are tabulated and plotted in Figure 6 with colors. The black 
lines in Figure 6 show the Z e -R relationships for falling snow presented in previous 
studies such as Sekhon and Srivastava (1970), Fujiyoshi et al. (1990), Vasiloff et al. 
(2000), Boucher and Weiler (1985), and Skofronick- Jackson et al. (2004). 


Comparisons in Figure 6 show that NWS radar observed Z e versus retrieved R 
curves of all snow crystal models considered in this study agree well with previously 
published Z e -R relationships. Retrieved snowfall rates assuming the C4 snow particle 
model, yield the smallest values, while retrieval results based upon HCs yield the largest 
snowfall rate for a given radar reflectivity. Variations of retrieved snowfall rates can 
differ by a factor of 2. Further, the ranges of uncertainties associated with different 
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particle shapes enclose all previously published Z e -R relationships and all snow particle 
shapes evaluated in this study. 

5.2 Consistency Checks 

The distribution of fractions of snow cover on the ground selected by the Bayesian 
retrieval method are shown in Figure 7(a). Within the AMSU-B field of view the snow 
cover fraction ranges between 0.2 and 0.6 over the snow storm core region. Figure 7(b) 
shows that the uncertainties (standard deviations) of the snow cover fraction are about 
0.15 to 0.2 in the storm core. These higher uncertainties are likely caused by the fact that 
none of the channels probe to the surface in the heavily precipitating storm core, and 
hence the Bayesian technique can select wide range of surface snow cover fraction values 
without affecting the resultant Tbs. 

It is noted that retrieved fractions of snow cover over the Great Lakes and St. 
Lawrence River regions (near 73°W and 46°N) are large. This is a direct result of the 89 
GHz channel sensitivity to surface features. The 89 GHz AMSU-B image 48 hours earlier 
than the time analyzed for the snowfall retrievals and prior to the start of our analyzed 
snow fall event is shown in Figure 8(a). In this image, a cold surface feature is seen in the 
same region where the retrieved snow cover fraction is high. Indeed, all the AMSU-B 89 
GHz images from March 1-10, 2001 exhibit the same cold brightness temperatures. This 
region is a river valley with elevations near sea level. To validate of retrieved surface 
snow coverage in this region, snow water equivalent (SWE) values were derived with the 
algorithm in Foster et al. (2005) using the Special Sensor Microwave Imager (SSM/I) 19 
GHz and 37 GHz vertically polarized channels. The SSM/I retrieved SWE for March 4, 
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2001 are shown in Figure 8(b), prior to the blizzard analyzed in this work. The SWE 
distribution provides evidence that there was more snow on the ground near the Great 
Lakes and St. Lawrence River regions than in other areas before the March 5-6, 2001 
blizzard, suggesting that the distribution of retrieved surface snow fraction derived in the 
current study is valid. 

The relative humidity (%) and uncertainties at 0.02 km altitude above the surface in 
the retrieved profiles are shown in Figure 9. In Figure 9(a), retrieved relative humidity 
values are large (above 90 %) over the snow storm core region where the radar 
reflectivity is greater than 30 dBZ (Fig. 1(a)). The uncertainties range up to 7 %. The 
uncertainties are lower over the region where the retrieved relative humidity values are 
high. It should be emphasized that the uncertainty here is in reference to the MM5 
produced relative humidity profiles. Therefore, the accuracy of the retrieved humidities 
is limited by the range of humidities in the profiles in the MM5 data base. 

Compared to the previous retrieval results of Skofronick- Jackson et al. (2004), the 
new algorithm retrieves snowfall rate distributions that are more consistent with the NWS 
radar reflectivity distribution (Fig. 6). It should be noted that the distributions of retrieved 
parameters such as snow water content and relative humidity in Skofronick-Jackson et al. 
(2004) were noisy, and that the maximum retrieved snow water content was associated 
with physically inconsistent low relative humidity profiles in that study. In the current 
study, the improved physical and statistical assumptions in the retrieval algorithm result 
in relative humidity values that are relatively high over the snowstorm region, which is 
consistent with moisture convergence and lifting in this intense storm (Fig. 9). 
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6. SUMMARY AND CONCLUSIONS 


In this study, we improved a snowfall retrieval method by addressing our 
assumptions described in Skofronick-Jackson et al. (2004) and estimate snowfall rate 
over land using the AMSU-B Tb observations at 89 GHz, 150 GHz, 183. 3±1 GHz, 
183.3±3 GHz, and 183.3±7 GHz. This improved retrieval model is applied to the 
blizzard of March 5, 2001 over New England. Major improvements resulting from this 
study are following: 

• The current study employs the Discrete -Dipole Approximation (DDA) 
method (Purcell and Pennypacker 1973; Draine and Flatau 2003) to calculate 
single scattering parameters for various nonspherical snow particles in 
radiative transfer calculations. 

• This study incorporates snow particle size distributions based on in-situ 
observations by Houze et al (1979), Lo and Passarelli (1982), and vary with 
height. 

• Instead of assuming that the fall velocity is fixed at 1 m/sec for all falling 
snow particles (Skofronick-Jackson et al. 2004) to convert retrieved snow 
water content to snowfall rate, this study employs the relationships between 
fall velocity and snow particle sizes described in Rutledge and Hobbs (1983). 

• The current study uses a Bayesian method for optimizing the retrievals and 
correlations between variables are considered in retrievals. Unlike previous 
Bayesian retrievals of Noh et al. (2006), the present study accounts for the 
error covariance of Tbs in retrievals. This reduces the noise that was present 
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• This study demonstrates that modeling errors are strongly correlated 
especially between 89 GHz and 150 GHz channels and between 183.3±3 and 
183. 3±1 GHz channels. Water vapor uncertainties result in Tb errors that are 
significantly anti-correlated at the lower and higher channel frequencies. The 
largest Tb error covariances are associated with the 89 and 150 GHz 
channels, primarily due to assumed uncertainties in modeled surface 
emissivities and PSDs. 

Comparisons of the current retrieved snowfall rates with NWS radar reflectivity 
measurements indicate better consistency in relation to those of Skofronick- Jackson et al. 
(2004). This consistency is evaluated in terms of reduced noise in the retrieved 
distributions, the use of non-spherical snowflake shapes, correlations between variables, a 
Bayesian inversion method with Tb error covariances between different channels, and the 
use of realistic surface emissivity variables. Results suggest that the physical model 
developed in this study uses more appropriate and realistic assumptions and further 
improves estimates of snowfall rate over land based upon high frequency microwave 
brightness temperatures. 

The results of error analysis show that the uncertainty inherent in the retrievals of 
snowfall ranges between 0.01 mm/hr and 2 mm/hr at 0.02 km altitude above the surface, 
increasing with snowfall rate. For retrieved relative humidity, the uncertainty is less than 
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7 % over the whole retrieval domain, although this uncertainty is more of a measure of 
the error in the variability of the MM5 relative humidity profiles. Standard deviations of 
retrieved surface snow coverage over most of the storm area range between 0 and 0.2. 


APPENDIX: Sensitivity 

The sensitivities of Tbs to variations of humidity, hydrometeor profiles, and surface 
emissivity, represented by Jacobians, helps to explain the model error covariance derived 
in section 4: Greater sensitivity of Tbs to a given geophysical parameter means that there 
is greater model error associated with uncertainties in that parameter. In addition, the 
analysis of model Jacobians is standard procedure in data assimilation schemes because 
Jacobians provide crucial information on model sensitivity to input perturbations as a 
function of model state. Its application to cloud and precipitation parameters via radiative 
transfer analysis is relatively new (Moreau et al. 2003). It has also influenced radiometer 
optimization studies outside the NWP community (Bauer and Mugnai 2003). Here, 
Jacobians associated with relative humidity and hydrometeor profiles, which are 
extracted out of the MM5 simulations (Section 3.1) over land, are calculated. 

Figure Al and Figure A2 show Jacobians, J, for AMSU-B channels with respect to 
humidity and snow water content profiles, respectively. The J units for humidity and 
snow water content profiles are AK/A(g/kg) and AK/A(g/m ), respectively. In these 
figures, solid and dashed lines indicate mean values and standard deviation values, 
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respectively. The sensitivity of the calculated Tbs to the fraction of snow coverage and 
snow particle size distributions is presented in Figure A3. The snow particle model used 
to calculate Tbs shown in these figures is HC (Fig. 2). 

From Figure Al, it may be noted that the 89 GHz channel is nearly insensitive to 
water vapor variations, and that the 150 GHz Tbs slightly increase with water vapor 
increment (< 0.5 K/(g/km)). As we may expect, the water vapor channels (183.3 ± 1, 
183.3 ± 3, and 183.3 ± 7 GHz) are sensitive to variation in humidity. In particular, the 
calculated Tbs are most sensitive to variations of water vapor amount at altitudes around 
7 km - 8 km. Considering the water vapor channels, on average, the most opaque channel 
(e.g., 183. 3±1 GHz) is more sensitive to water vapor variations than the least opaque 
channel (e.g., 183.3±7 GHz). 

In Figure A2, the calculated Jacobians show that both the 183.3 ±7 GHz and 150 
GHz channels are most sensitive to snow water content. It is also noted that the water 
vapor channels are not sensitive to the snow water content near the surface, where the 
calculated Jacobians at 89 GHz and 150 GHz show sensitivities of 1 K/(0.1 g/nr) and 2.5 
K/(0.1 g/m 3 ), respectively. The 89 GHz, 150 GHz, 183.3 ± 1 GHz, and 183.3 ± 3 GHz, 
and 183.3 ± 7 GHz channels are most sensitive to the snow water content at 3 km, 4 km, 
5 km, 6 km, and 7 km altitudes, respectively. This result suggests a limited snow 
profiling capability if these channels are combined in a retrieval algorithm and the water 
vapor connection can be determined or specified independently. 
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Figure A3 presents histograms of Tb differences (ATb) when the fraction of surface 


snow coverage increases from f = 0.5 to f = 0.6. From the figure it is seen that the 89 
GHz channel is three times more sensitive to the surface snow coverage than the 150 
GHz channel, while the water vapor channels are nearly insensitive to surface emissivity 
changes. 
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Figure 1. (a) Radar reflectivity (dBZ) obtained from the NWS operational radar composite at 
variable heights between 0.5 and 2.5 km measured at 2300 UTC on March 5, 2001. Note that 
heavy snowfall over CT, MA, NH, and VT. Tbs (K) observed from (b) 89 GHz, (c) 150 GHz, (d) 
183.3 ± 7 GHz, (e) 183.3 ± 3 GHz, and (f) 183.3 ± 1 GHz of AMSU-B at 23:02 UTC on March 
2001 . The color scale shown in (b) applies to the subsequent figures. 
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(a) (b) 




Figure 3. Profiles of (a) wind speed and (b) Richardson number measured in Brookhaven, NY at 
18 UTC on March 5 and 00 UTC on March 6, 2001 during the New England blizzard. 
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Figure 4. Retrieved (a) water equivalent snowfall rate (mm/hr) and (b) uncertainty estimate 
(mm/hr) at 0.02 km altitude. Since the retrievals are applied to precipitation over land, the oceanic 
1 regions are masked. 
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Figure 5. Departures of estimated brightness temperatures from measured brightness 
temperatures: Tb(estimated) - Tb(measured) (K) at (a) 89 GH,z (b) 150 GHz, (c) 183.3 ± 7 GHz, 
(d) 183.3 ± 3 GHz, and (e) 183.3 ± 1 GHz. 
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Figure 6. Measured NWS radar reflectivity Z e versus retrieved snowfall rate (mm/hr) using 
different snow particle models. Results are compared with previously published Z e -R relationships 
for falling snow. The uncertainty values of snowfall rates retrieved different snow particle models 
are shown as error bars. 
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Figure 7. Retrieved (a) fraction of snow cover on the ground and (b) its uncertainty estimate. 
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Figure 8. (a) The 89 GHz AMSU-B channel image on March 3, 2001. (b) Distribution of snow 
water equivalent (SWE) derived with the algorithm in Foster et al. (2005) using the SSM/I 19 GHz 
and 37 GHz vertically polarized channels, over the area of interest on March 4, 2001 . 
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Figure 9. Retrieved (a) relative humidity (%) and (b) uncertainty estimate at 0.02 km altitude. 
Since the retrievals are applied to precipitation over land, the oceanic regions are masked. 
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Figure A1 . Jacobians for AMSU-B channels with respect to water vapor mixing ratio. 
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Figure A2. Jacobians for AMSU-B channels with respect to snow water content (swc). 
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Figure A3. Histograms of calculated Tb differences when the fraction of surface snow coverage 
increases from f = 0.5 to f = 0.6. 
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